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Abstract 

This work is a continuation of our previous investigation of binary data 
corruption due to a Brownian agent [T. J. Newman and W. Triampo, preprint 



cond-mat/9811237 ]. We extend our study in three main directions which 
allow us to make closer contact with real bistable systems. These are i) a 
detailed analysis of two dimensions, ii) the case of competing agents, and 
iii) the cases of asymmetric and quenched random couplings. Most of our 
results are obtained by extending our original phenomenological model, and 
are supported by extensive numerical simulations. 
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I. INTRODUCTION 



This paper is the second part of a two-stage investigation into the statistics of an active 
random walker (Brownian agent) in a bistable medium. This is but one example of the 
myriad of systems in which a random walker interacts in some way with its environment 
We consider systems in which the Brownian agent (BA) performs a pure, unbiased 
random walk in a medium composed of elements which may take one of two values. On 
visiting a given element the BA has a certain probability of switching the value of that 
element. In our first paper B (hereafter referred to as DCI) we motivated our investigation 
into such processes using the example of data corruption caused by a Brownian agent (i.e. 
the elements were taken to be bits of binary data) and we outlined the possible applications 
of our results to describing soft error production in small-scale memory devices HH. The 
disturbance of a bistable medium by a BA is also related to reversible chemical kinetics by 
a high mobility catalyst, and disordering of bi-atomic structures by a wandering agent 0,0, 
such as an anion or cation vacancy in NaCl, or an impurity in a semiconductor compound 
(e.g. Zn in GaAs). One may also view this process as a non-conserved cousin of magnetic 
disordering via spin exchange with a wandering vacancy [|3|-|T0|j. Last, but not least, the 
analysis of this process yields a deeper understanding of the statistics of random walks. 

In DCI we studied the simplest possible process, namely a single BA disordering a bistable 
system, with a switching probability which is independent of the value of the element. One 
of our main concerns in DCI was to construct a phenomenological continuum model of the 
process. This enabled us to find results independent of microscopic details, and also to 
study coarse-grained quantities, such as the probability distribution for the local "density 
of disorder" which is less easily defined on a lattice. Our results were derived mostly for the 
case of one spatial dimension. 

In order to make contact with a wider range of processes it is necessary to consider a more 
general model. This is the aim of the present work. We shall extend our original investiga- 
tion in three directions. First, we shall present a careful analysis of two spatial dimensions. 
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We find that the phenomenological model predicts the correct asymptotic behavior, despite 
the need to introduce some form of regularization. We present results for the mean density of 
disorder p, and using some special properties of the continuum description, we shall also de- 
rive an approximate form for the probability distribution of the density of disorder. Second, 
we shall consider a system containing more than one BA, thereby inducing "competition" 
as each BA interferes with the disorder created by the others. The main result here is that 
the disordering efficacy (i.e. the global amount of disorder due to N agents as compared 
to one agent) is massively reduced for dimensions less than two, whereas in precisely two 
dimensions, each BA eventually becomes independent (in that the disorder it creates is not 
reordered by other BA's). We shall present calculations based on the continuum model to 
make these statements quantitative. Third, we shall consider two kinds of generalized cou- 
plings between the BA and the bistable medium: asymmetric switching probabilities, and 
quenched random switching probabilities. We shall argue that these generalized couplings 
may be modeled within the continuum limit by simple generalizations of our original model, 
and we shall derive some basic consequences, which for the case of quenched randomness 
are particularly interesting (see the outline below). In all cases, we shall support our results 
by numerical simulations of the underlying lattice model. 

A more detailed outline follows: In section II we present a recapitulation of the results of 
DCI. We shall briefly describe the lattice model in general spatial dimension d, and present 
it in the form of a master equation. We then describe the associated continuum theory, and 
we state without derivation some pertinent results previously derived from this continuum 
model in DCI. In section III we concentrate on calculating p in two spatial dimensions. Given 
that two is the critical dimension of the process, the results are modulated by logarithmic 
corrections. Therefore great care is needed to compare different theoretical predictions and 
numerical results. We shall present (briefly) four alternative methods of calculation, and 
show that they all predict the same asymptotic behavior, and agree with the numerical 
simulations once the strong corrections to scaling are included. Using these results, we are 
also able to approximately reconstruct the probability distribution for the density of disorder. 
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In section IV we study the case of more than one BA in the context of a generalized version 
of our continuum model. We shall derive an exact integral expression (for asymptotically 
large times) for the disordering efficacy cr N (d) for arbitrary N. First we concentrate on 
d — 1. We evaluate ctjv(1) for N = 2,3,4, and also extract its functional form for N 3> 1. 
We then extend our study to arbitrary spatial dimension d < 2, and evaluate 0" 2 (<f) and 
(fN^d) for large N. We use these results to analytically continue to two dimensions, thereby 
avoiding the use of a microscopic regularization. In section V we consider the system in 
d — 1 with generalized couplings. First we study asymmetric rates, so, for example, in the 
data corruption process the BA will have different probabilities to switch — ► 1 and 1 — > 0. 
We propose to model this using a simple extension of the original continuum theory, based 
on the idea that relative to the non-zero "background disorder" the dynamics of the system 
are the same as the symmetric case. Second, we consider quenched random couplings. In 
this case we argue that p picks up logarithmic corrections in time, while the global amount of 
disorder remains unaffected. We shall also discuss the quenched average of the distribution 
function of disorder density, and show that it is very sensitive to the distribution of the 
couplings. We shall support our results by numerical simulations which are described in 
detail in section VI. We end the paper with a summary of our results and our conclusions. 



II. RECAPITULATION 

In this section we give a very brief review of the main ideas and some of the results 
contained in DCI in order to place the present work in a proper context. The process of a 
BA in a bistable medium is first modeled on a hypercubic lattice of dimension d. The position 
of the BA is denoted by a lattice vector R. In a time step St the BA has a probability p to 
move to one of its 2d nearest neighbor sites. In making such a jump, there is a probability 
q that the element on the site departed from is switched. The elements are described by 
spin variables a r (where r denotes a discrete lattice vector) which may take the values ±1. 
The spin variables encode the information about the disordering process. For example in the 
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data corruption process we label uncorrupted bits (of value 1) by spin +1 and corrupted bits 
(of value 0) by —1. [Thus, we shall often use the terms "magnetization density" and "global 
magnetization" which may be simply translated as "density of disorder" and "total amount 
of disorder".] This process is illustrated in Fig.l for d = 2 and p — q — 1. We can define 
the dynamics via the probability distribution P(R, {a r },t), which is the probability that at 
time t, the BA is at position R and the spins have values given by {cx r }. This distribution 
evolves according to a master equation Jll| which takes the form 



P(R, K}, t + 5t) = (l- p)P(R, {a r }, t) + P{l 2d q) £ P(R + 1, K}, t) 

+ |j£P(R + l,---,-<TR + i,---,t) , (1) 

where {1} represents the 2d orthogonal lattice vectors (which have magnitude I). 

The most direct quantities to extract from the master equation are marginal averages, 
the simplest of which is the magnetization density given by 

e(n, R, t) = Tr ff a ri P(R, R}, t) . (2) 

Averaging the master equation over the spin variables gives 

0(r, R, t + St)- 6(r, R, t) = ^ £[6(r, R + 1, t) - 6(r, R, t)} 

-^e(r,r,0EW- ( 3 ) 

At this stage a continuum limit (in both space and time) may be taken of the above 
equation, which yields 

$9(r, R, t) = ^V^6(r, R, t) - A 0(r, R, t) A,(r - R) . (4) 

Two parameters have appeared: the effective diffusion constant of the BA given by D = 
2l 2 p/5t, and an effective coupling between the BA and the spins given by A oc pql d /5t. 
This continuum equation for has the form of a diffusion equation with a sink potential 
Ai(r) which is a strongly localized function with lateral extent / and normalized to unity. 
In the naive continuum limit, this function may be taken to be a d- dimensional Dirac delta 



function. However, for d > 2 it is necessary to smear this function in order to regularize the 
theory. 

In DCI an alternative continuum description was obtained by viewing the process as a 
stochastic cellular automaton. The process is then defined in terms of the position R(£) of 
the BA (which is now an independent stochastic process), and the coarse-grained density of 
disorder (or magnetization density) which is defined in a small region of space at a specific 
time, and is a functional of R(t). In some sense, one may view this in the same spirit 
as a Langevin description of a stochastic process described at a more fundamental level 
by a master equation. Taking the continuum limit of this description yields first a simple 
Langevin equation for the position of the BA 

<^R j— y \ ,-\ 

■5- = £(«), (5) 

where £(£) is a noise term, each component of which is an uncorrelated Gaussian random 
variable with zero mean (i.e. ^(t) is a white noise process). The correlator of £ is given by 

Mt)tiV)) = B'6 i j6(t-lf). (6) 

Here and henceforth, angled brackets indicate an average over the noise (or equivalently the 
paths of the BA). The BA is chosen to reside initially at the origin: R(0) = 0. 
The evolution of the magnetization density is described by 

d t <f>(r, t) = -A>(r, t) A t (r - R(t)) . (7) 

This equation may be integrated to give the explicit functional solution 



(f)(r, t) = exp 



-X' J dt' A,(r-R(0) 
o 



(8) 



The above solution is obtained for an initial condition 0(r, 0) = 1, which we shall use 
exclusively. In terms of the original lattice model it corresponds to choosing all the spins 
to have the initial value of +1, so that we measure the subsequent disorder of the system 
by counting the number of minus spins in the system. [Although of no relevance to the 
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coarse-grained description, we mention here that it is often convenient to choose the spin at 
the origin (i.e. the initial BA position) to be —1, so that for p = q = 1 all spins have value 
+ 1 after the first jump of the BA.] 

In DCI we showed that the continuum descriptions (|) and (|7|) are indeed equivalent 
with the identifications D = D' and A = A'. [This was proven by considering the former as 
an imaginary time Schrodinger equation Hl2| , and writing the solution of the latter as an 
imaginary time Feynman path integral [131-1 This ends our review of the modelization of 



the process - the reader is referred to DCI for further details and discussions. 

In DCI we exclusively used the continuum description (H) to generate results for various 
average quantities. The simplest quantity to consider is the mean magnetization density 
given by m(r, t) = (0(r, t)) (which is equivalent to the sum of G(r, R, t) over the BA position 
R). For d < 2 there is no necessity for regularization and we replaced the sink function A; 
by the d- dimensional Dirac delta function. The (temporal) Laplace transform for m(r, t) 
was found to have the exact form 



m(r, s) = - 

s 



\g{r,s) 



(9) 



1 + A£(0,s) 

where g(r, s) is the Laplace transform of the diffusion equation Green function 

g(r,t) = (2nDt)- d/2 exp(-r 2 /2Dt) . (10) 

This exact result allows one to extract a great deal of statistical information about the 
process. First, one may simply invert the Laplace transform to find the average magnetiza- 
tion density (or average density of disorder relative to 1/2) as a function of r and t. Explicit 
forms are given in DCI for d — 1. We note here that the average magnetization density at 
the origin decays for long times in d — 1 as 



2D \ 1/2 

m(0,t) 



1 + 



D 
A^ 



(11) 



,7rA 2 t, 

The continuum solution has the important property that (0(r,t;A) n ) = (<f)(r,t;n\)). This 
allows us to utilize the exact solution (^) to reconstruct the probability density for the 
magnetization density. We define V via 



P(0,r,t) = (5(0-0 R (r,t))) , (12) 

where 0R,(r, t) is the stochastic field solution given in Eq. (||). As explained in DCI, one 
may solve for this distribution exactly. In particular, for d = 1 the probability distribution 
for the magnetization density at the origin takes the form 



W,0,i) = ^iexp 



(log(0))< 



(13) 



4XH 

which is a log-normal distribution (and where we have defined A = \/ (2D) 1 / 2 ). This is 
interesting, as it indicates the extreme nature of the fluctuations in this system. For instance, 
the typical value of the magnetization density can be found from the above expression to 
decay exponentially, i.e. typ ~ exp(— A t/2), whereas the mean density decays as 1/Av^ as 
given in ([TT]) . 

Another interesting quantity which may be extracted from m(0, t) is the average global 
magnetization defined (relative to its initial value) as 

M(t) = Jd d r [<0(r,O)>-(0(r,t)>] . (14) 

As shown in DCI, for d < 2 this quantity obeys the exact relation 

^ = M<M). (15) 

Thus, for large times in one dimension we have M{t) ~ (Dt) 1 ^ 2 independent of the coupling 
A. In other words, the total amount of disorder created by a single BA on average increases 
as (-Dt) 1 / 2 , and is (rather surprisingly) independent of the coupling between the BA and the 
spins (for times larger than D/X 2 ). 



III. TWO DIMENSIONAL SYSTEMS 

In this section we shall present a careful analysis of the case of two dimensions. The 
simple random walk is recurrent for dimensions d < 2, whereas for d > 2 the walker has a 
probability less than unity for ever returning to its starting point Q. This basic fact from the 
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theory of random walks has an immediate implication for our data corruption problem. The 
non-recurrent nature of random walks for d > 2 implies that the BA will continually corrupt 
new regions of the system, and rarely revisit sites which it has previously corrupted. Thus 
the relative density of disorder (or average magnetization density) at the origin m(0, t) will 
decay to a non-zero (and non-universal) value, and the total amount of disorder (or average 
global magnetization) M(t) will increase linearly in time, with a non-universal prefactor. 
For practical applications, in which one wishes to limit the disordering capabilities of the 
BA, the first requirement is to restrict the geometry of the system to a dimension d < 2. So 
d — 2 is the critical dimension of the problem, and because of this we can expect logarithmic 
corrections to modulate the leading order results, and also to cause long cross-over times, 
thus making numerical results more difficult to interpret. 

In DCI we studied some general properties of higher dimensional systems and we also 
derived an approximate form for m(0, t) in d = 2 using a crude form of regularization. 
Rough agreement was found between this result and simulations, but no quantitative data 
analysis was performed. In this section we shall rederive the form of m(0, t) more carefully. 
The reasons for this are threefold: first, so that we can have confidence in the leading order 
result, and also get some idea of the sub-leading corrections; second, to provide insight into 
the relation between the discrete and continuum approaches in two dimensions (which is 
important given that the latter must be regularized for d — 2); and third, to allow us to 
construct the form of the probability distribution of the magnetization density (for which 
we need as much information about m(0,t) as possible). 

We shall use four different methods (which will each be described with brevity) to derive 
the form of m(0,t). Each has its strong and weak points as we shall see. The methods 
to be used are i) calculation from the exact lattice formulation for the marginal average 
(0), ii) solution of the diffusion equation (|J) using a smeared sink function, iii) calculation 
with infinite-order perturbation theory of the continuum theory (|8]) using a crude temporal 
cut-off (as described in DCI), and finally iv) analytic continuation from the exact result ([|) 
valid for d < 2. 
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A. Lattice calculation from Eq.(j3|) 

Referring to the equation of motion (|3]) for the marginal average we set d = 2, and for 
convenience we set the hopping probability p — 1, giving 

9(r, R, t + 5t) = \ ]T 9(r, R + 1, t) - | 0(r, r, t) £ 5 r , R+I . (16) 
4 l ^ l 

As an initial condition we take the BA to be located at the origin, and all spins to be 

+ 1, except the spin at the origin which is taken to be —1. [This convention is useful for 

q = 1 so that all spins have value +1 after one time step. However if q <C 1 it is more 

convenient to take the spin at the origin to be initially +1 as the chance of it flipping 

after one time step is small. We stress that these different choices for the initial value of 

the spin at the origin have no effect on the asymptotic properties of the system, and only 

serve to smoothen the magnetization density in the immediate vicinity of the origin.] Thus 

9(r, R, 0) = <5r,o(1 — 2<5 rj0 ). The solution of fll6|) may be attained by discrete Fourier and 

Laplace transform. Defining the former via 

e(r,M = £e(r,R,t)e ik - R (17) 

R 

and the latter via 

oo 

0(r, R,z) = J2 ( r > R > n5t) * n , ( 18 ) 

?1=0 

it is fairly straightforward to diagonalize Eq.(|T6"D to the form 

A fr k z) _ (l-26r, )-2qzf(k)e**Q(T,T,z) 

U(r, k, z) - 1 _ s/( k ) ' ^ iJ J 

where /(k) = [cos k x l + cos /c 2 /]/2. One may now solve the above equation self-consistently 
for 0(r, r, z) by inverting the discrete Fourier transform. One has 

V[t,t,z) 1 + 2qz Jdkf(k)D(k,z) ' [Z[)) 

where -D(k) = [1 — zfik.)}^ 1 and the momentum integrals are over the two-dimensional Bril- 
louin zone. The average magnetization density at the origin is given by summing O(0, R, t) 
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over R, which is equivalent to the zero Fourier mode 0(0, 0, t). Thus, substituting (p0| ) into 
(pT9|) we have after some rearrangement 

~l + q'(l -z)Jdk D{k,z) 



X)e(o,R,z) 



1 + q' Jdk D{k,z) 

where q' = 2g/ (1 — 2g). [The function J dk -D(k, z) is very well known in the theory of ran- 
dom walks 0] , and is the discrete Laplace transform of the probability of a random walker 
to return to its starting point after n steps.] Finally we must inverse Laplace transform the 
above equation. For large n we can extract the asymptotic form of the average magnetiza- 
tion by invoking the Tauberian theorem In this case we need the form of the Laplace 
transform as z — *■ 1. Using |I[ 

J dk D(k, *) ~ ~ log (j^j [1 + 0(1 - z)\ , (22) 

we have from the Tauberian theorem the asymptotic result 

Y 6(0, R, nSt) = , , . ~ l . (23) 

^ (qt'/tt) log n + C(q) 

The constant C(q) is not accessible from the Tauberian theorem, although it could in prin- 
ciple be calculated from a careful inverse Laplace transform of Eq . (f21p . 

In our simulations we generally take the switching probability q to be unity which implies 
q' = — 2 and thus 

g e <°- = e/,) icgn- car (24) 

In section VI we shall make a direct comparison of this result with numerical simulations 
in two dimensions. It is also interesting to note from Eq. (|2i|) that setting the switching 



probability q = 1/2 gives 2©(0, R, ndt) = —S n0 . In other words, the average value of the 

R 

spin at the origin remains at zero after one jump of the agent. This "maximal uncertainty" 
for q = 1/2 only holds at the origin, and spins at other lattice sites will have a positive mean 
for all times. As a final note, if q 1 it is more convenient to choose the initial value of the 
spin at the origin to be +1 in which case we find from the foregoing analysis 

y 6(0, R, nSt) = j- . - 1 — — . (25) 

r (2g/vr) log n + C{q) 
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B. Diffusion equation (^) with smeared sink function 



We now wish to solve for the coarse-grained magnetization density m(0, t) within some 
continuum limit. In this subsection we accomplish this by solving the diffusion equation 
(|j) with a smeared sink function. This calculation is the closest in spirit to the lattice 
calculation since (^) is the direct continuum analog of the discrete equation (|^) solved in 
the previous subsection. Setting r = 0, Eq.(fJ) takes the form 

«9 t 0(O, R, t) = ^V^0(O, R, t) - A 6(0, R, f)A,(R) . (26) 

The simplest finite range form to take for the sink function is a radial function which is 
zero outside a radius /, and (through normalization) equal to (l/irl 2 ) within this radius. It 
is also convenient to smear the initial condition of 6 in a similar way: 



( 1/tt/ 2 , 


|R 


< I 


l». 


|R| 


> I 



0(O,R,O) = (27) 



Given the radial symmetry of both the sink function and the initial condition, 0(0, R, i) 
is a function only of the radial coordinate R and t. Thus we define Oi(-R, t) and &2{R,t) 
to describe the original function respectively inside and outside the radius /. On Laplace 
transforming these functions in time, it is straightforward to show using boundary value 
techniques that 

8) = nl2 ^ +XI) + A ^ (V 7 ^ R') 

& 2 (R,s)=B(s)K ( v ^R') , (28) 



where we have introduced the scaled quantities R' = Ri^l/D) 1 ^ 2 and A' = \/nl 2 , and where 
Io{z) and Kq(z) are modified Bessel functions fi~4j] . 



The functions A(s) and B(s) can be determined by demanding continuity of and its 
radial derivative at the boundary R = I. Their forms are complicated and we shall not write 
them explicitly here. The average magnetization density m(0, t) is obtained by integrating 
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0(0, R, t) over R. Using the properties of integrals of modified Bessel functions [14| we have 
the closed form expression 

/ x 1 2A'C0) , > 

m(0 -'> = ITv + W»(/+v)v» - (29) 

where 



c(s) = wn/i^n 

and V = 1(2/ D) 1 / 2 . We may now take the lateral size I of the sink function to zero, and 
retain the leading terms in the above expression. Using the properties of Bessel functions 



at small argument we finally arrive at 

1 



m(0, s) = 

s 



(31) 



_w(X) - X\og(sl' 2 ) 

where A = X/2nD and w(x) = y / x/o(2\/a;)//i(2y / x) which approaches unity for small x. It 
is important to mention that this expression is not valid for Re (s) > l/l', and thus the 
apparent pole is an artifact of the limit I — ■> 0. Examination of the complete form ( |55| ) 
indicates there is no pole for Re (s) > 0, as expected on physical grounds. We refer the 
reader to Appendix A in which this Laplace transform is inverted. For large t the result 
takes the form 

m(0,t) = — B 1 , (32) 

w{X)+X\og(t/r) 

where r = / /2 e -7 and 7 = 0.57721... is Euler's constant. 

Thus we have arrived at our result, and indeed the leading order behavior has the same 
functional form as Eq.(|25|) derived in the previous subsection. A closer comparison of the 
leading terms indicates that A oc q for q 1 as expected on physical grounds. In the present 
continuum calculation, we have also extracted an explicit form for the sub-leading term (i.e. 
w(X)) which is non- universal and depends on the precise form of the smeared sink function. 
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C. Infinite order perturbation theory of Eq. (|8|) 



In DCI we concentrated on analyzing the continuum solution Eq. (|^) using infinite order 



perturbation theory in A, using a Dirac delta function in place of the sink function Aj. For 
d < 2 each term may be averaged, and the series may then be resummed exactly after 



in two dimensions a similar procedure may be used, but now each term in the perturbation 
expansion diverges after averaging over the noise. A controlled analysis requires one to work 
with a particular sink function, but as noted in DCI, resummation of the perturbation series 
is very difficult in this case. Instead, we persisted in using the Dirac delta function as a sink 
function, and regularized each term in the perturbation theory by introducing a short-time 
cut-off to, the physical origin of which is the microscopic correlation time of the noise. This 
regularization scheme is not systematic, and this is one of the reasons for the present more 
careful analysis. We shall enter into no details here of the perturbation theory method as a 
comprehensive description may be found in DCI. We also refer the reader to Ref. || where 
the same regularization scheme was used in calculations on vacancy-mediated diffusion, and 
produced scaling results in accord with exact lattice calculations ||. 

Using the time cut-off regularization scheme allows one to extract the dominant contri- 
bution from each term in the perturbation expansion, and to resum the series. The result 
is 



It is interesting to compare this result with Eq . (|32"D obtained in the previous subsection. 
They are seen to agree (for small A, in which case w(X) ~ 1) if we make the identification 
to = T {= 2/ 2 e~ 7 / D). This indeed supports the role of to as a microscopic correlation time 
of the noise, since it is seen to correspond to the time taken for the diffusion process to 
correlate a microscopic region of size ~ I 2 . 



Laplace transform. No divergences appear, and one arrives at the exact result (^). Exactly 




(33) 
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D. Analytic continuation from Eq.@ 



As a fourth method of extracting the form of m(0, t) for d = 2 we briefly mention analytic 
continuation from d < 2. As discussed in the previous section, no regularization is required 
in the perturbation expansion method for d < 2, and one retrieves the exact result 

1 



m(0, s) 



1 + Xg(0,s) 



(34) 



where 



9(0, s) 



r(l -d/2) 
(2nD) d / 2 s l - d / 2 



and V{z) is the gamma function [T4[. Using the result r(^) ~ l/z for 2 
e = 2 - d -> 



m(0, s) 



1 



(35) 

we have for 
(36) 



.1 - Alog(s/s ). 
where Sq = e _2 / e . 

On comparing this result with Eq. (|3l|) obtained in section IIIB, we see that they are 
equivalent (for small A, in which case w(X) ~ 1) if we make the identification 1/e = log(l//'). 
This relation has only a formal meaning, as there is no physical sense in continuing dimen- 
sionality. However, it is instructive to learn that the leading order result may be retrieved 
intact from analytic continuation. Indeed, we shall use this tool in section IV to study 
competing agents in two dimensions. 



E. Calculation of the probability distribution 

To complete this section, we shall briefly discuss the probability distribution of the 
magnetization density at the origin. As mentioned in section II, the continuum theory (|5p 
has the property (0(0, t; A) n ) = (0(0, t;n\)). Thus, knowledge of the A dependence of the 
first moment allows one to reconstruct all the moments of the magnetization density, and 
thus the probability density V((j),0,t) for this quantity. The method for retrieving V from 
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the moments is given in detail in DCI for the case of d = 1. A similar calculation suffices 
for d — 2, so long as one has an "accurate" form for the mean magnetization density. In the 
present section we have attacked this problem from four different directions and have arrived 
at agreement for the asymptotically dominant term for m(0, t). However, as noted in DCI, 
it is necessary to have more information than this in order to construct V correctly. If one 
tries to calculate V using m(0, t;X) ~ 1/Alogt one obtains a distribution V ~ m(0,t)/(f) 
which is singular at = 0. The sub-leading correction to m(0, t; A) is crucial to determine 
the distribution correctly. Our results from the previous four subsections are in agreement 
with regard to the time dependence, but differ in the way in which A appears in the sub- 
leading term. We prefer Eq.(|32|) in that an explicit form w(X) appears. However, we have 
been unable so far to reconstruct V using this form due to the complicated nature of w. For 
small A, w(\) ~ 1, and Eq.(|32|) then coincides with the less controlled results (|33|) and ( flop 
of subsections IIIC and HID respectively. The reconstruction of V is possible from these 
forms and one finds 

P( ( p,0,t)=P(t)^- 1 , (37) 

with (3{t) = l/\\og(t/to). As t — > oo this distribution approaches the form m(O,t)/0, but 
is never singular for finite times. In section VI we shall describe our attempts to measure 
V for two dimensional systems. Our numerical results are in surprisingly good agreement 
with Eq. (|37|) above. 

IV. COMPETING AGENTS 

In this section we shall analyze the effects of many BA's within the system. We shall 
assume the BA's to be non-interacting, in the sense that they are unaware of each other's 
immediate presence. The non-trivial statistics reside in the fact that the disordering effects 
of the BA's statistically interact via the overlap of the BA histories. As we have already 
seen, a single BA interferes with the previous disorder it has created, such that the amount 
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of disorder does not simply increase linearly in time. This effect is more severe when more 
than one BA is present, as each BA can disturb the disorder that another BA has previously 
created. 

We measure the strength of this interference via a quantity called the "disordering effi- 
cacy" of the agents, defined as 

MW(t) 



(38) 



where M N (i) is the average global disorder created by N agents. If the BA's were truly 
independent (in terms of the disorder they create) then we would expect <7n = N. As we 
shall see, for d < 2 the value of a N is strongly reduced below this value. However, for d = 2 
this value is recovered, but only in the deep asymptotic regime (t ^> e N ). 

The extension to many BA's is easily modeled within the continuum theory. We introduce 
N random walkers, each of which is described by a position vector TL a (t), a = 1, 2, • • • , N. 
Since the BA's are independent we have 



dt 



(39) 



where are independent Gaussian white noise sources with zero mean. The equation of 
motion for the coarse-grained magnetization density 0^ takes the form 

N 



d t ^ N \v,t) = -X^ N \r,t) £ A,(r - R a ( t )) , 



(40) 



a=l 



with solution 



A? 



0^(r,t)= []exp 



z 

-A J dt' A^r-R^t')) 



(41) 



On averaging over the paths of the iV agents, we have the particularly simple result 



m (N \r,t) = (0 (A °M)> = m (1) (r,t) 



A? 



(42) 



Thus the average global magnetization is given by 

M {N \r,t) ee J d d r [mW(r,0)-mW(r,t)] = J d d r 



l-m^(r,t) N 



(43) 
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Scaling the spatial coordinate by the diffusion length scale (2Dt) d ^ 2 , we have from Eqs. ( |38|) 
and © 

Jd d r \l-m^(r/(2Dt) 1 / 2 ,t) N ] 
a N {d) = lim — - L M „ or ^ 1/c , * ■ (44) 



t-oo fd d r [1 -m( 1 )(r/(2 J Dt) 1 /2 ) t)] " 
To proceed with the calculation it is convenient to first perform the large- 1 limit. We 
concentrate on d < 2, and therefore replace A/(r) by the Dirac delta-function. From Eqs. 
(§) and (|10D it is a fairly straightforward matter to show that 

T 

2 

t-+oo T(l ~d/2) 

We may use this result to explicitly evaluate the denominator of Eq. ([44]) giving us 



m(r) = lim m (1) (r/(2£>£) 1/2 , t) = — — - / du u l ~ d e"" 2 . (45) 

*^oo Tl 1 — a/2) J 







<7jv(d) = d T(l - d/2) J dr r d - x (1 - m{r) N ) . (46) 

o 

We refer the reader to Appendix B where it is demonstrated that the above expression may 
be recast in the more useful form 



DC 



aN ( d) = 2N j N _jj / dr r i-d e -*- m (r) N - 2 , (47) 



o 



r(l-d/2) 

for N > 1. We note that the result for two agents follows immediately from this expression, 
and we have o<i = 2 d l 2 . This result is striking. For the case of d = 1, we see that two 
agents create only \/2 as much disorder as one agent. Also, assuming we may continue this 
result to exactly two dimensions, we find that <7 2 (2) — 2 - i.e. two agents create disorder 
independently. In the remainder of this section we shall use Eq.([|7]) for two purposes. First, 
we shall concentrate on d = 1 and evaluate the exact values of cr 3 (l) and cr 4 (l), which may 
be used to compare with numerical simulations. Second, we shall evaluate the integral for 
large for arbitrary d G [0, 2] using a saddle point method. This calculation will make 
clear the tremendous difference in the large- A behavior of crjv(d) for d < 2 and d = 2. We 
end the section with a simple scaling argument which helps us to understand these analytic 
results. 

For d = 1 the expression fl4"5] ) is simply m(x) = erf(x), where erf(z) is the error function 



14j| . Thus we have from Eq. 



Ml) = 2N %2 ^ / dx e" 2 ^ erf(x)^" 2 



7T 



(48) 



We are able to evaluate these integrals for N = 3 and N = 4. The details may be found in 
Appendix C. The results are 



(7 2 (1) = V2 , ^(1) = ^ sin" 1 



(7 4 (1) = sin - 

7T \0 



(49) 



The numerical values of these expressions are presented in Table I along with the results 
from our computer simulations (the details of which may be found in section VI) . Excellent 
agreement is found. 



N 


Cjv(l)theory 


Cjv(l)simul 


2 


1.414- • • 


1.42(1) 


3 


1.662- • • 


1.66(1) 


4 


1.835 - - - 


1.84(2) 



Table 1: The predicted values of ctjv(1) for N = 2,3 and 4, from Eq . (f49|) , compared with 
our numerical simulations. 

We now return to the case of arbitrary d G [0,2] and consider the limit of large N. For 
the sake of generality, we study the integral 



Q N (f3,d) = J dr r x ~ d e' Pr m{r) N , 
o 

with m(r) as given above in Eq. (f45|) . We can recover the disordering efficacy via 

. , 2N(N-1) . . 

We wish to implement a saddle point calculation, so we rewrite (pOl) as 

Q N (P,d) = J dr r l ~ d e~ FN{T) , 

Q 

where 



(50) 



(51) 



(52) 



F N (r;P,d) = (3r 2 - N\og(m(r)) 



(53) 
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The saddle point is denned via dF^ / dr\ r=ro 
equation 



0, which yields for tq the transcendental 



& - Ne ' Tl 

1 m(r )r(l - d/2) ' 



(54) 



For N 1, it is easy to see that r ~ (logiV) 1 / 2 ^> 1, so that F(r ) ^> 1 and the saddle 
point method is self-consistently justified. There are two stages to the calculation. First we 
must evaluate F(r ) to the desired precision. Second, we must Taylor expand F(r) about 
the saddle point to account for fluctuations. [It turns out that we need to evaluate the 
leading and sub-leading terms in order to have a useful comparison to the numerical data. 
Thus the calculation is somewhat involved.] For the first stage we iteratively solve Eq.(|54"l) 
to obtain 



log 



N 



d 2 log(logiV) 
4 log N 



1 



logiV 



J3T (d/2) (log N) d / 2 
Substituting this result into Eq . fl53|) , we have after some manipulations 
~/3r(l-d/2)(logiV)< i/2n/3 



(55) 



exp[-F N (r Q )] 



eN 



1 ^ 2 log(logiV) | Q 
4 logA^ 



1 



logA^ 



(56) 



The second stage proceeds as follows. We Taylor expand the function Fn about its value 
at the saddle point: 

(r - r ) n d n F(r) 



F N (r) = F N (r ) + ^ 



n=2 



dr r 



(57) 



r=rg 



When using the saddle point method, it is usually only necessary to consider the second 
derivative in the above expansion (the so-called Gaussian fluctuations). In the present case, 
it turns out that all terms in the expansion contribute equally to the fluctuations. It is 
relatively straightforward to show by repeated differentiation of F^(r) that 



d n F(r) 



dr r 



P{-2r Q 



r=ro 



1+0 4 



(58) 



Setting r = r + f we have from Eqs. (|57| ) and ( ^ ' 

-2r f) T 



F N (r) 



F N (r )+PYl 



n=2 



ni 



O 



(59) 
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We now substitute this expansion into Eq.([5^). Scaling f by 2ro and explicitly performing 
the sum over n we have 



Q N ((3,d) = — L / dr exp 



2r$ 



-(3 (e~ f + f - l) 



-2ri 



(60) 



Performing the integral and neglecting exponentially small terms we have 



1 + Ok 



(61) 



Finally, we substitute into this expression the form fl55|) for r , and the result ( |56"D for -F7v(r ) 
to yield 



Qn{M 



Y(f3)Y{l-d/2f {\o g N)^ d ' 2 
2 



(/?-l)rf 2 log(logiV) 
4 log 



1 



logiV 



(62) 



Combining Eqs. (|51~1) and (^) then gives us the final large- N result 



a N (d) = T(l - d/2) (log N) 



d/2 



] rf 2 log(logAQ | Q 



1 



(63) 



4 logAf \logN y 
In particular, for d = 1 the disordering efficacy increases as (log N) 1 ^ 2 with strong logarithmic 
corrections. We also note that for real- valued d < 2, the disordering efficacy increases as 
(logiV) d//2 . In other words, the BA's overlap very strongly in their disordering for all d < 2. 

The main reason for calculating <y^{d) for real- valued d < 2 was to attempt to analytically 
continue the result to d = 2 which is otherwise difficult to calculate due to the presence of an 
explicit regulator. However, if we simply set d = 2 in our expression for a^(d) above, we find 
that the apparent behavior 0^(2) ~ logiV is multiplied by the infinite constant. This gives 
us the hint that the large- N behavior for d = 2 is stronger than logiV, but offers no more 
information than that. We can, however, trace the failure of the saddle point calculation 
for d = 2 back to the saddle point equation ( p>4] ) for r . For large- iV and d < 2 we have the 
leading result ro ~ [logeiV] 1 / 2 , where e = 2 — d (and we have obtained this result from ( |54]) 
by taking the small e and large r limit of m(r)). Thus for fixed e > we may always take 
iV large enough to create a saddle-point at a large value of r , in which case the result ( |63|) 
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is valid. However for fixed iV (no matter how large), taking e — > squeezes the saddle point 
back into the origin, in which case the saddle point method is of no use. 

Therefore we cannot use Eq.(^) to analytically continue to d — 2. However, a simpler 
method may be used to extract the result. We refer the reader to the original expressions 
for a^{d) and m(r) as given in Eqs. Q4"7|) and Q4"5| ) respectively. Taking e C 1 and N ^> 1 
these equations take the form 

oo 

a N (d) = eN 2 J dr r e ~ l e~ 2r2 m(r) N , (64) 
o 

and 

r 

m {r) = ejdu u e ' x e~ u2 . (65) 
o 

On studying Eq.flBip we see that as e — > the integral is dominated by small r. Thus 
we require the small-r form for m(r) which is easily extracted from ( p5| ) to be m(r) ~ r € . 
We now break the integral in (|4]) into two pieces. The first piece encompasses the range 
r G (0, 1) so that we can neglect the Gaussian factor and substitute the small r form of m(r) 
to find the leading order result 

i 

cr^(d) = e N 2 J dr r Nt ~ l = N . (66) 
o 

The second piece encompasses the range r 6 (1, oo) and for large enough N will have the 
asymptotic form cr^(c?) ~ (log AT) /e. To summarize, the contribution crj^(d) dominates for 
C 1/e, and the contribution ^(d) dominates for N ^> 1/e. Thus for d = 2 we have the 
result <7at(2) = iV for all N. This is consistent with the analytic continuation to d = 2 of 
the exact result cr 2 (d) = 2 d l 2 found earlier. 

These results may be understood in the following way. Consider iV agents in a system of 
d dimensions. Since each agent performs a random walk, the amount of available space in 
which we can expect to find the agents has a volume ~ t d l 2 . Also, we know from our previous 
results that the amount of disorder created by a single agent increases as t d l 2 for d < 2, 
but as t/logt for d = 2. Thus for d < 2 the available space for A^ agents and the amount 
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of disorder created by a single agent both scale as t d l 2 which means that there is certain 
to be interference between the agents, which will become severe for large iV as we have 
seen. However for d = 2 the amount of available space for N agents scales logarithmically 
faster with t than does the amount of disorder created by a single agent. Thus for times 
large enough such that N <C logt we expect the disorder created by the agents to become 
statistically independent and thus (7^(2) = N. Note that the "independence time" grows 
exponentially with the number of agents, making numerical observation difficult for even 
moderate values of N. Similar arguments have been made in the context of "the number of 



distinct sites visited by N random walkers" |16|]. In section VI we shall present results for 
iV = 2, 3 and 4 in two dimensions. The disordering efficacy is seen to slowly approach its 
expected value. 



V. GENERALIZED COUPLINGS 

In all of our work so far we have taken the coupling between the BA and the spins to be 
symmetric and spatially homogeneous. At least one of these properties is likely to be absent 
in a practical application of our model. Asymmetry in the switching probability is likely 
in data corruption since the states of a bit are not physically encoded in a symmetric way. 
Also in chemical kinetic applications, the reaction rates between two chemical species are 
unlikely to be symmetric. Spatial homogeneity of the couplings is also an idealized situation. 
Generally there will be quenched random fluctuations in the switching rates of different sites, 
and it is important to quantify the effect of this randomness on the results obtained so far. 
In the next two subsections we shall consider these important generalizations in turn. 



A. Asymmetric rates 

Consider the lattice model described in section II with the additional property that the 
probability of flipping a spin depends on the value of the spin. If the BA leaves a site with 
spin +1 we flip that spin with probability q + , whereas if the BA leaves a site with spin — 1 
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we flip that spin with probability q~ . At a microscopic level the model is now considerably 
more complicated as the transition rates for flipping depend explicitly on the values of the 
spins. We shall not construct the master equation for this case. Rather, we shall try to 
construct an analog of the continuum model (0) using a simple physical idea. 

The main effect of the asymmetric couplings is to favor one type of spin over the other. 
For the sake of argument let us take q + < q~ . In this case the probability that a site with 
spin +1 is flipped by the BA is less than that for a site with spin —1, so that after some 
time there will be considerably more sites in the active zone with spin +1 than sites with 
spin —1. By "active zone" we mean the region which has been thoroughly explored by the 
BA after an elapsed time t. This region will have a linear size which grows as \f~Dt. We may 
immediately obtain a rough estimate of the average global magnetization from this picture. 
Let us make the crude approximation that outside this region the magnetization density 
relative to its initial value of unity is zero, and that inside this region the relative density is 
1 — m cq . The quantity m eq is the equilibrium value of the average spin value which has the 
form 

m eq = _ . 67 

q + q + 

We may therefore estimate that the average total magnetization with asymmetric rates is 
related to that with symmetric rates via 

M A {t) = (1 - m ccl )M s (t) = {1 + ^ /q+) M s (t) , (68) 

where the subscripts A and S indicate "asymmetric" and "symmetric" respectively. [Using 
this simple argument for the case q + = q~ would imply that M{t) ~ (Dt) 1 ^ 2 , independent 
of the rates since the quantity m eq = 0. This is indeed the described in section II. 1 

We now make the approximation that in a coarse-grained model the dynamics of relax- 
ation to this "background magnetization" are the same as the relaxation to zero magnetiza- 
tion in the symmetric model. This amounts to a simple linear shift in the order parameter 
4>(r, t) which appears in the continuum description. We therefore write the analog of (0) for 
asymmetric rates in the form 
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d t <f>{r, t) = -A[0(r, t) - m eq ] A,(r - R(t)) . (69) 

We may now use this model to derive results for the average magnetization density and 
the average global magnetization. Define the field 5<p(r,t) = <p(r,t) — m cq . Then 5<p(r,t) 
satisfies the original continuum theory (0), but with an initial condition 1 — m eq . Thus, for 
any realization of the BA we have the relationship 

(f) A (r, t; 1) - m eq = 5 (r, t; 1 - m eq ) , (70) 

where we have included the initial value of the field as the final argument. Given that the 
field is linearly proportional to its initial value, the above relation simplifies further to 

4> A (r, t; 1) = m cq - (1 - m eq )0,s(r, t] 1) . (71) 

Thus the average magnetization density relaxes to its equilibrium value with the dynamics 
of the symmetric system (cf. Eq. (|TTD ) , up to a factor of (1 — m eq ). Integrating this equation 
over space, and using the definition (|1^) of the average global magnetization, it is easy to 
see that the relation ( |B"BD is an exact consequence of this model. 

We have tested the validity of ([71]) by measuring the average magnetization density of a 
patch of spins in a numerical simulation, and have found reasonable agreement. We report 
on our numerical work in more detail in the next section. 

B. Quenched random rates 

We now consider an alternative generalization in which the rates are symmetric, but 
spatially inhomogeneous. At the microscopic level, this is modeled by attaching to each 
lattice position r a quenched random variable q r G (0, 1) drawn from some distribution 
S({q r }). The random variable q r gives the probability of switching in the event that the 
BA visits the site r. In the continuum theory, we model these quenched random couplings 
by a simple extension of Eq.(|l]). We generalize the coupling parameter A to a spatially 
inhomogeneous random function A(r), with a distribution S[X\. This does not affect the 
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solvability of the model due to the locality of the continuum theory, so that we have the 
explicit solution 



(r, t; A) = exp 



-A(r) J dt ! A,(r-R(0) 
o 



(72) 



where we have emphasized the dependence of <fi on A. 

There is an enormous variety in the forms of the quenched disorder distribution S that 
one may study. A comprehensive analysis is beyond the scope of the present work. In our 
simulations (to be described in the next section) we have chosen to study switching rates q r G 
(0, 1) which are independent random variables drawn from a uniform distribution. There is 
more than one way in which such a distribution can manifest itself at the continuum level. 
Depending on the extent to which one coarse-grains the lattice model, the continuum theory 
will have a distribution of couplings which is either very similar to the lattice distribution (i. e. 
uncorrelated and uniform), or more Gaussian in nature (due to the central limit theorem). 
We expect different physics according to the behavior of the distribution S[X] as A — > 0. 
Roughly speaking, if S has zero weight for small- A we expect the quenched disorder to be 
irrelevant to the system for large times. Conversely, if S has non-zero weight for small A, 
then the disorder will play a role for arbitrarily late times. In order to exemplify this latter 
case, we consider as an example the distribution 

S[A] = nSioc(A(r)), (73) 

r 

where S\ oc is a local (or "on-site") distribution function which we take to be uniform: 

5i oc (A) = (1/A m ) H(X m - X)H(X) , (74) 



where H(x) is the Heaviside step function ||17|| . Henceforth we shall concentrate on d = 1 
for simplicity. 

Due to the way in which the quenched random couplings enter the continuum theory, we 
can easily generalize our earlier exact results. Let us first concentrate on the magnetization 
density at the origin, averaged over both the BA trajectories, and the distribution of A. The 
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first average is performed in the usual way using Laplace transforms (see DCI) and we have 
the analog of Eq. @ 

1 



m(0, s; A) = - 

s 



1 + X(0)g(0,s)_ 

On averaging this density over the distribution of random couplings (|74]) we have 



(75) 



/ 2D \ 1/2 
(m(0,s;\)) s = ^ log 



1 + 



Ai\ 1/2 ' 



(76) 



2Ds J 

Inverse Laplace transforming yields the asymptotic result 

< m(o - (;A M^r iog (^) • (7?) 

Therefore the decay in the homogeneous case (of 1/y/t behavior) for the average magnetiza- 
tion density is slowed by a logarithmic factor due to the presence of the quenched random 
couplings. It is interesting that this logarithmic slowing down does not affect the average 
of the global magnetization. This quantity may be explicitly calculated using the method 
of infinite order perturbation theory described in DCI. The result is that the leading order 
term for large times is independent of the coupling, and is unaffected by the average over S. 
Thus the average over the random couplings will in no way affect the leading order behavior 
of M{t) ~ (Dt) 1 / 2 . The average of the global magnetization is a very robust quantity, and is 
unaffected by scalings of the homogeneous couplings (as seen in section II) and by making 
the couplings quenched random variables. We have performed numerical simulations for 
systems with quenched random couplings and report our findings toward the end of the 
following section. 

Finally, we briefly discuss the changes that can occur to the probability distribution 
V of the local magnetization density at the origin, when quenched random couplings are 
introduced. We have from section II the exact result ( |T3D for this function in d — 1: it is a 
log-normal distribution. The same form will hold in the present case, but now the parameter 
A is to be averaged over S[X\. We average over the uncorrelated uniform distribution given 
above in Eq.(fT~|). Simple integration gives 
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(V)s = * WflSM*!), (78) 

where A m = X m / (2D) 1 ^ 2 and Ei(z) is the exponential integral |L4] . This distribution behaves 
very much like the pure log- normal distribution for <C 1. However, we have the interesting 
result that the distribution is singular at <fi = 1. Explicitly we have for ip = 1 — <fi -C 1: 

2(x()V2A ra 

This logarithmic singularity stems from the fact that if S[X] has a non-zero weight for 
arbitrarily small A, there will be a non- vanishing subset of systems for which almost no 
corruption occurs. For large times, this subset will appear as a pronounced "peak" in the 
distribution (V)s for — > 1. 



VI. NUMERICAL SIMULATIONS 

In this section we describe the details of our numerical simulations, and present figures 
of our data to support the various theoretical claims made in the preceding sections. The 
simulations are performed on a rf-dimensional lattice (with d either 1 or 2) and follow the 
lattice rules described in section II. We set the hopping probability p = 1, so that in each 
time step the BA is moved randomly to one of its 2d nearest neighbor sites. In doing so, 
the spin it leaves behind is flipped with probability q (which is set to unity unless otherwise 
stated). Systems are taken large enough so that the BA never touches the boundary during 
the lifetime of the simulation. We generally average over between 10 5 and 10 7 realizations, 
depending on the desired precision of the simulation. 

In section III we calculated the asymptotic form of the average magnetization density at 
the origin from an exact lattice calculation. As noted, the sub-leading term is not obtainable 
from the Tauberian theorem. However, it is clear from the various calculations in section III 
that the corrections decay as l/(logt) 2 . We have checked this in the following way. From our 
simulations we measure m(0, £) ( = J2n ©(0, R, n8t) ). We take its inverse and subtract the 
predicted asymptotic result of (2/tt) logn (in order to compare with Eq.fl24|)). The resulting 
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data is plotted against 1/logn, which should yield a straight line, the ^-intercept of which 
fixes C(l). As seen in the inset of Fig. 2, the data indeed confirms this expectation, and we 
find C(l) ~ —0.742. In Fig. 2, we have plotted the inverse of the data against logt, along 
with the predicted asymptotic result. The large difference between the curves indicates the 
strong role of the corrections (which are of order 10% even for t ~ 10 6 ). 

We have also attempted to measure the probability distribution V(<f), 0, t) in d = 2, as 
calculated in section HIE. In DCI we attempted to measure this quantity in d = 1 and we 
met with limited success. The difficulty lies in the fact that the distribution only makes 
sense for a coarse-grained magnetization density (since a single spin always has a bimodal 
distribution). Thus, to construct V numerically we must measure the magnetization density 
for a patch of spins. If the patch is too small, the resulting histogram will have too few bins 
to have any smooth structure; but if the patch is too large, the time scales required to allow 
the agent to wander far away from the patch and return many times become numerically 
prohibitive. Thus we must compromise, and we use a square patch of 121 spins. The agent 
is started on the boundary of the patch to avoid the patch being internally decimated by 
the transient motion of the agent. The predicted result for V is given in Eq. (|3"T|) . It shows 
three types of behavior set by a cross-over time t* defined by (3{t*) = 1. For t <C t* the 
distribution has most weight near = 1, and then for t ~ t* the distribution becomes 
uniform over the entire interval of 0. For late times t ^> f, the distribution approaches 
the form ~ 1/0. Our numerical measurement of V is shown in Fig. 3 for these three time 
regimes. Good qualitative agreement is found. [Note that the measured distribution has 
non-zero weight for negative due to the modest patch size, and is also suppressed near 
= 1 due to internal decimation of the patch magnetization.] 

In section IV we presented an analysis of the continuum theory with N agents. We 
described the interference between the agents by a number termed the "disordering efficacy" 
(JN{d). In particular we calculated the exact value for o"at(1) for N = 2, 3 and 4, as given 
in Eq. (f4"9"D . We also predicted the large- N form of a^(l), along with the exact form of 
the strong logarithmic corrections, as shown in Eq.(^). In two dimensions we argued that 
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<7jv(2) = N for all N, but only for times t 3> e . We have performed numerical simulations 
of the many agent system in order to test these results. The microscopic rule we use is that 
there is no hard-core exclusion between the agents, and that for each time step the N agents 
are in turn moved to a randomly chosen nearest neighbor site. A spin which is occupied by 
two agents, say, will thus (for q — 1) be flipped twice in that time step. In Fig. 4 we show the 
evolution of the ratio of the average global magnetization for N agents as compared to one 
agent, for d — 1. Asymptotically this ratio is the disordering efficacy by definition. Results 
are shown for N = 2, 3 and 4. The curves asymptote to constants as expected, the values 
of which are compared to the theoretical predictions. Excellent agreement (to within less 
than 2%) is found. The numerical values are also given in Table 1 in section IV. We have 
measured the disordering efficacy for higher values of N in d = 1, and we plot our results on 
a logarithmic scale in Fig. 5. Also shown is the theoretical prediction (^) for large- N (where 
we have included both the leading and sub-leading terms). Again, excellent agreement is 
found. As argued in section IV, the cross-over times in the two-dimensional many agent 
system are very large, growing exponentially with N. We have attempted to measure <Jn{2) 
numerically, but we have been unable to reach large enough times to see the ratio of the 
average global magnetizations reach a constant value. However, we may still compare our 
results to the theoretical prediction, by plotting the predicted value of the asymptote (= N) 
minus the measured ratio, against 1/logt. If the predicted asymptote is correct, the data 
so plotted should converge to the origin. In Fig. 6 we present such a plot for N = 2, 3 and 4. 
The resulting curves appear to be heading for the origin, and therefore constitute numerical 
support for the prediction a^{2) = N. 

In the first part of section V we considered a system with asymmetric switching rates, 
and argued for the very simple phenomenological description given by the equation of motion 
(|69|). This equation may be solved and one obtains the simple relationship ( fTTD between 
the magnetization density for asymmetric rates and that for symmetric rates. In particular, 
this result may be used to derive the intuitive result ( p%| ) for the proportionality of the 
average global magnetization for the asymmetric and symmetric cases (Ma(£) and Ms(t) 

30 



respectively). We have performed simulations of the system with asymmetric rates in d = 1. 
In Fig. 7 we show the ratio of M^(t) /Mg(t) for different choices of the rates q + and q~. It 
is seen that the ratio gradually tends to a constant, and that the constant is in agreement 
with the theoretical prediction (|68|). We have also tried to test Eq. ([7TD at the level of the 
average magnetization density. We have found that the ratio of the densities is constant 
for all but the shortest times, which already indicates that the form of the time decay 
of m(0,t) (relative to its equilibrium value) is insensitive to the asymmetry in the rates. 
However, we have found that this constant is sensitive to whether we measure the densities 
for an individual spin, or for a patch. The value of the constant decreases as we increase 
the patch size. In Fig. 8 we plot the ratio of magnetization densities for asymmetric and 
symmetric systems, for various patch sizes with the rates q + = 0.8 and q~ = 0.2. As we 
see, the ratio seems to be approximately constant, and the value decreases as the patch size 
is increased. Reasonable agreement is found with the prediction of Eq.([71]) for the largest 
patch of 101 spins. This indicates that ([H]) is likely to be correct, but only after a large 
degree of coarse-graining. 

Finally we turn to the predictions of section VB, which concern the system with sym- 
metric quenched random rates. We investigated this situation numerically in d = 1 by 
measuring the global magnetization, and the magnetization density. We perform the aver- 
aging in batches; namely, we use the same set of quenched rates for a batch of 1000 systems, 
which are then averaged over their different BA histories. Then we repeat this for Nf, batches 
(with N b ~ 10 3 ), thus averaging over different realizations of the quenched rates. We choose 
the quenched couplings {q r } to be uncorrelated random variables drawn from a uniform 
distribution in the range [0, 1]. In Fig. 9 we show the average of the global magnetization in 
a system with quenched rates, and compare it to the same quantity in a "pure" system. The 
two curves become identical for late times (following a power law growth ~ indicating 
that the average global magnetization is insensitive to the quenched rates, as predicted. 
Turning to the average magnetization density, we found at first that our results could not be 
fitted to the theoretical prediction (|77|) . However, a closer analysis revealed that great care 
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must be taken in making the comparison. The point is that one wishes to compare theory 



and simulation for long times (when the asymptotic form given in (|77 ) becomes valid). How- 
ever, it turns out that there is a new cross-over time tf within the numerical simulations, 
beyond which the theoretical prediction is expected to break down. This cross-over time 
emerges due to insufficient averaging over batches. In the simulation we average over Nf, 
batches, and thus select Nj, couplings from the uniform distribution. There will be a smallest 
selected value qf ~ l/iVj, of the coupling, and thus from this finite sampling of S, we cannot 
discern if we are really using a uniform distribution with non-zero weight all the way down 
to q = 0, or a distribution with non-zero weight only down to q — qf. In the latter case, 
one can shown (from either the lattice or continuum theories) that the logarithmic slowing 
down will vanish after a time tf ~ Thus, in our simulations we can only expect to 

see the logarithmic slowing down for times much less than tf ~ N£. So, in order to makes 
this time-window as large as possible it is important to perform as many batch averages 
as possible, at the expense of averaging over BA histories within a given batch. To make 
a quantitative comparison to theory, we have calculated the exact form of the logarithmic 
slowing down from the lattice theory, using the methods described in section IIIA. After 
averaging over a uniform distribution of the rates, we have from the Tauberian theorem the 
exact asymptotic result 



O 



(80) 



log 71^ 

We plot in Fig. 10 the average magnetization density, and the above asymptotic result. 
Good agreement is found, and the difference may be well fitted against the slow logarithmic 
corrections. Note, for times longer than those shown on the plot, the data fails to agree with 
the predicted result due to cross-over into the regime t ^> tf in which finite sampling effects 
dominate. 
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VII. SUMMARY AND CONCLUSIONS 



This paper has been devoted to a continuation of our previous study DCI |3| of the 
statistical properties of data corruption due to a Brownian agent (BA). In order to make 
closer contact to potential physical realizations of this process we have extended our original 
study in three directions: i) a careful examination of two dimensions, which is the critical 
dimension for the model, ii) the competition induced by more than one agent, and iii) the 
case of generalized couplings between the BA and its environment. 

In section II we provided a recapitulation of DCI and outlined the discrete and continuum 
modelizations of the process, along with presenting results pertinent to the present study. 

Section III was devoted to a careful examination of two dimensions, which, because of the 
recurrent properties of random walks, is the critical dimension for this process. Thus, naive 
dimensional arguments are modified by logarithmic corrections, and there are also strong 
corrections to the leading order results. For these reasons it is difficult to obtain controlled 
analytic results and to compare them meaningfully to numerical data. In an attempt to 
overcome these obstacles we presented four different methods of calculation for the average 
magnetization density at the origin m(0,t), using both discrete and continuum models. We 
found agreement in all cases for the leading behavior m(0,t) ~ 2nD/\logt. We also gained 
insight into the non-universal nature of the sub-leading terms, which in all cases depends 
on the type of regularization used. However it is clear from the various forms of m(0,t) 
that the leading corrections decay as l/(logt) 2 . This allowed us to analyze our data in 
a quantitative way, and excellent agreement was found with our numerical simulations for 
both the leading term, and for the nature of the logarithmic corrections. Furthermore, given 
the explicit functional forms derived for m(0,t) we were able to derive approximate forms 
for the probability distribution V of the magnetization density at the origin. The precise 
form is again non-universal as the shape of V depends sensitively on the sub-leading terms 
of m(0,t). However, we compared the simplest form fl37|) with simulation and found good 
qualitative agreement. 
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In section IV we generalized the model to a process with N agents, which although 
non-interacting, produce statistical interference in the corruption of the system due to the 
overlap of their histories. We described their effect via the "disordering efficacy" 0jv((i) 
defined as the ratio of the average global magnetization due to N agents, as compared to 
that of one agent. After deriving a general form for a N (d) (for d < 2), we first studied 
the case d — 1. We were able to derive exact values for <7jv(1) for N = 2,3 and 4 which 
are in excellent agreement with simulations. Returning to general d < 2 we found that 
02 (d) = 2 d / 2 , which assuming smooth continuation to two dimensions implies 02(2) = 2, 
meaning that the two agents are statistically independent. We then formulated a large- N 
analysis of the disordering efficacy, and found for d < 2 that a N (d) ~ [log N] d / 2 . Thus 
in dimensions less than two, the agents interfere very strongly with one another. The 
precise logarithmic dependence (with the accompanying strong logarithmic corrections) are 
supported by excellent quantitative agreement with simulations in d — 1. We found that 
the large- N analysis fails on continuation to d — 2. An alternative analysis for d — 2 
yielded the result 0tv(2) = N which is satisfied in the deep asymptotic regime t^>e N . This 
result is also supported by numerical simulation in d — 2 for modest values of N . Thus, 
the disordering efficacy is extremely sensitive to the dimensionality of the system, and is 
strongly discontinuous as d /* 2. 

We moved on to a study of generalized couplings in section V. First we considered 
homogeneous couplings (between the BA and the environment) which are asymmetric. In 
the context of data corruption, this means that the agent has a probability of switching 
— > 1 which is different to the probability of switching 1 — ► 0. We proposed to model this 
effect by a simple extension of the original continuum model, based on the idea that the 
dynamics of the asymmetric process may be described relative to the non-zero background 
magnetization, using the relaxational dynamics of the symmetric process. This idea leads 
to several simple predictions for the average (local and global) magnetization which are 
supported by simulation. Second, we considered symmetric but inhomogeneous couplings. 
Thus, the probability for the agent to switch a given element is given by a quenched random 
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variable. We incorporated this new effect into the continuum model by generalizing the 
constant parameter A to a quenched random field A(r) described by a distribution S[X\. We 
investigated the simple case in which 5* is spatially uncorrelated, with a uniform on-site 
distribution function. We found that the average magnetization density at the origin decays 
more slowly than in the homogeneous case, by a factor ~ log t. However, the average global 
magnetization is asymptotically insensitive to the quenched couplings. The probability 
density for the magnetization density at the origin V(<f),0,t) is dramatically altered for 
values of close to unity, where it suffers a logarithmic divergence; this is due to a non- 
vanishing subset of systems for which the coupling X(x = 0) is very small, such that barely 
any corruption occurs. The main features of corruption in the presence of quenched random 
couplings were supported by numerical simulation. 

Finally in section VI we described the details of our numerical simulations, and presented 
figures supporting the different analytic predictions from the previous sections. 

In conclusion, we have presented three important extensions to our original work. The 
most important fact to emerge from this investigation is the robustness of the phenomeno- 
logical model (^) in accounting for the properties of the data corruption process under a 
wide range of conditions. The fact that this process may be modeled by such a simple con- 
tinuum theory gives one confidence in applying similar phenomenological models to more 
complicated physical situations involving Brownian agents. 

As mentioned in the Introduction, the process of disordering of a bistable medium by a 
Brownian agent has putative applications to data corruption, homogeneous catalysis and the 
effect of impurities in biatomic crystals. It remains to be seen whether one can find a solid 
application of the models studied in DCI and the present work. By having explored more 
challenging situations such as two dimensions, many agents, and asymmetric and quenched 
random couplings, we are able to seriously address the practical issues involved in making the 
connection between the rich statistical properties of our model, and their potential existence 
in the real world. 
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APPENDIX A: 



In this appendix we outline the procedure of taking the inverse Laplace transform of 
Eq.(|3TD. Consider the function H(t), the Laplace transform of which has the form 

1 



His) 



(Al) 



s [a — b log s] ' 

for sufficiently small Re(s) such that the apparent pole at s — e a / b may be disregarded 
(a and b are positive constants). The inverse Laplace transform is given by the Bromwich 
integral 



H(t) 



ds 



,,st 



c 2iri s [a — b log s] 



(A2) 



where the contour C is to be taken up the imaginary axis (note the integration around the 
simple pole at the origin gives a contribution of zero). The important singularity is the 
branch point at the origin. We take the cut along the negative real axis and deform the 
contour around this cut. Integrating across the cut then gives us 



H(t) 



dx 



b J x 
o 



-xt 



(A3) 



(logx — a/b) 2 + 7r 2 

We are interested in the asymptotic form of this integral for large t. In this limit the integral 
is dominated by small x. We therefore expand the denominator to give 



H(t) 



dx 



-xt 



bJ^ x (logx) 2 



2a 

1 + 7^ + 



1 



(A4) 



61ogx ' \^(logx) 2 

By splitting the range of integration into two pieces (x G (0, 1/t) and x G (l/£, oo) ) we 
systematically evaluate the terms in the above expansion to find 



H(t) 



(a + 76) 



+ 



(A5) 



blogt (Mogt) 2 ' ~ \{logt) 3 J ' 

where 7 = 0.57721... is Euler's constant. Using the appropriate form for the constants a 
and b we arrive at Eq.(|32"D correct up to 0(l/(logt) 3 ). 
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APPENDIX B: 



In this appendix we indicate the steps leading from Eq.fljq) to Eq. (|4y|) . First we factorize 
the integrand to give 

N-1 

a N {d) = dF(l - d/2) £ Jn{d) , (Bl) 

?1=0 

where 



oo 



J n {d) = I dr r d ~ l f{r)m{r) n , (B2) 



with f{r) = 1 — m(r). We then integrate the above expression by parts using 

du = - + C0 " Sta " t • (B3) 

2 

The resulting integral contains a factor of re~ r which allows one to integrate by parts a 
second time, thus yielding after some algebra (for n > 0) 

Ud) = ~ dT(l-d/2) + Kn{d) ~ Kn ~ l{d) ' (B4) 



with 



/ , oo 

2n{n + 1) f x _ d 2r 2 ,. n _! 



K ^ = m^iw { irr e ™M • < B5 > 

Substitution of Eq.(|B4T) into (|Bl"l ) yields the result fl47D as shown in the main text. 

APPENDIX C: 

In this appendix we indicate the evaluation of 0jv(l) as given by Eq.(|48|) for N = 3 and 
N = 4. For iV = 3 we require 

oo 

<7 3 (1) = ^ j dx ^ erf (x) • ( C1 ) 


We introduce the integral 
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J(/3) = vr 1 / 2 | dx erf (x) . (C2) 

o 

Using integration by parts it is straightforward to show that 1(0) satisfies the differential 
equation 

2 4 + 1= -(Tffl • (C3) 

This differential equation may then be integrated (using the boundary condition /(oo) = 0) 
to give 

W = ^72 / V pfiw) = F? sin " ((TtW) ' (C4) 

Using this result along with the fact that cx 3 (l) = (12/7r)/(2), we have the leftmost result of 
Eq.(f49"D as shown in the main text. 

The evaluation of cr 4 (l) proceeds along similar grounds. We define a function 

oo 

L(0) = vr 1 / 2 J dx e-^ 2 erf (x) 2 , (C5) 
o 

which may be shown to satisfy the differential equation 

dL 2 
2 ^ + £ = - (l+fl(2 + flV (C6> 

This equation may be integrated, and the resulting integral may be calculated by elementary 
methods to yield 

L <« = ^ s " rl (rb) • (C7 > 

Using this result along with the fact that cr 4 (l) = (24/7r)L(2), we have the rightmost result 
of Eq. (f49"D as shown in the main text. 

A similar method may be used to simplify the integrals for higher N, but we have been 
unable to evaluate these integrals in closed form. Naturally, the integral in Eq.(^) may be 
evaluated numerically for a given iV to any desired precision. 
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FIG. 1. Illustration of the data corruption process for d = 2, with p = q = 1. The initial 
uncorrupted state is shown on the left, with the BA represented by the filled circle. On the right 
we show a typical walk of ~ 20 steps. The BA switches a bit with each visit, so those bits visited 
an even number of times are restored to their original value. 
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FIG. 2. Plot of l/m(0, t) (the inverse of the average magnetization density at the origin), as a 
function of \ogt for d = 2. The lower solid line is the asymptotic prediction, while the upper solid 
line (which is partially obscured by the data) is the prediction with corrections to scaling included. 
The inset shows the difference between l/m(0,t) and its asymptotic form, plotted against 1/logi. 
The data may be fitted to a straight line, thereby yielding the corrections to scaling explicitly. 
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FIG. 3. Plot of the probability distribution V as a function of cj) for d = 2, for times t = 10 3 
(diamonds), t = 10 4 (plusses), and t = 10 5 (squares). Note the qualitative similarity (as time 
proceeds) between these different histograms, and the theoretical prediction (| 
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FIG. 4. Plot of the ratio of the average global magnetization for iV agents as compared to a 
single agent, for (in ascending order) N = 2,3, and 4, in d = 1. The solid lines are the theoretical 
predictions (f49|) for the disordering efficacies, which are the asymptotes of this ratio. 
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FIG. 5. Plot of the disordering efficacy on as a function of N in d = 1. The solid line is the 
theoretical prediction (^) (including the strong logarithmic corrections). 
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FIG. 6. Plot of the ratio of the average global magnetization for N agents as compared to a 
single agent, subtracted from N, as a function of 1/logt for (in ascending order) N = 2, 3, and 4, 
in d = 2. The data in each case is apparently heading for the origin, thus supporting the theoretical 
prediction = N . 
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FIG. 7. Plot of the ratio M^(i) /Ms{t) versus time for systems in d = 1 with asymmetric rates. 
In ascending order, the values of the rates are (q~,q + ) = (0.5,0.2), (0.5,0.4), (0.4,0.5), (0.2,0.5). 
The solid lines are the theoretical predictions obtained from (|68|). 
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FIG. 8. Plot of the ratio of magnetization densities mA(0, t) /ms(0, t) for d = 1 with asymmetric 
rates q~ = 0.2 and q + = 0.8. The ratio is measured for different patch sizes, the patches having 
from top to bottom 9, 19, and 101 spins respectively. The solid line is the theoretical prediction 
from Eq.(0). 
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FIG. 9. Plot of the average global magnetization versus time for systems with quenched and 
pure (homogeneous) switching rates. For long times the functions become identical, growing ~ \/t. 
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FIG. 10. Plot of the magnetization density at the origin (averaged over both BA histories and 
the quenched rates) versus time. The solid line is the asymptotic form of the theoretical prediction 
Eq.(pO|) with no free parameters. 
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